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Abstract 


1 Introduction 


The Morse-Smale complex is an important tool for 
global topological analysis in various problems of 
computational geometry and topology. Algorithms 
for Morse-Smale complexes have been presented in 
case of piecewise linear manifolds [11]. However, 
previous research in this field is incomplete in the 
case of smooth functions. In the current paper we 
address the following question: Given an arbitrar¬ 
ily complex Morse-Smale system on a planar do¬ 
main, is it possible to compute its certified (topolog¬ 
ically correct) Morse-Smale complex? Towards this, 
we develop an algorithm using interval arithmetic 
to compute certified critical points and separatrices 
forming the Morse-Smale complexes of smooth func¬ 
tions on bounded planar domain. Our algorithm can 
also compute geometrically close Morse-Smale com¬ 
plexes. 


keyword Morse-Smale Complex, Certified Compu¬ 
tation, Interval Arithmetic. 


Geometrical shapes occurring in the real world are 
often extremely complex. To analyze them, one as¬ 
sociates a sufficiently smooth scalar field with the 
shape, e.g., a density function or a function interpo¬ 
lating gray values. Using this function, topological 
and geometrical information about the shape may 
be extracted, e.g., by computing its Morse-Smale 
complex. The cells of this complex are maximal 
connected sets consisting of orthogonal trajectories 
of the contour lines—curves of steepest ascent—with 
the same critical point of the function as origin and 
the same critical point as destination. The leftmost 
plots in Figures 1(a) and 1(b) illustrate the level sets 
of such a density function h, and the rightmost pic¬ 
tures the Morse-Smale complex of h. as computed 
by the algorithm in this paper. This complex re¬ 
veals the global topology of the shape. Recently, the 
Morse-Smale complex has been successfully applied 
in different areas like molecular shape analysis, image 
analysis, data and vector field simplification, visual¬ 
ization and detection of voids and clusters in galaxy 
distributions [6, 13]. 
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1.1 Problem statement 

A Morse function h. : —> R is a real-valued func¬ 
tion with non-degenerate critical points (i.e., crit¬ 
ical points with non-singular Hessian matrix). Non¬ 
degenerate critical points are isolated, and are either 
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(a) Contour plot (left) and Morse-Smale complex (right) 
of K(x,y) = cosx siny + 0.2 (x + y) inside box [—3.5, 3.5] x 
[—3.5, 3.5]. CPU-time: 11 seconds. 


(b) Contour plot (left) and Morse-Smale Complex (right) 
of h(x,y) = lOx — ^(x^ + y^) -|- (x^ -1- y^)^ inside box 

[—5,5] X [—5,5]. CPU-time: 0.5 seconds. 

Figure 1: Contour plots of Morse-Smale functions, 
and their Morse-Smale complexes. 

maxima, or minima, or saddle points. They corre¬ 
spond to singular points of the gradient vector field 
Vh. of h, of type sink, source or saddle, respectively. 
Regular integral curves of the gradient vector field 
Vh. are orthogonal trajectories of the regular level 
curves of h. We are interested in the configuration 
of integral curves of the gradient vector field. An 
unstable (stable) separatrix of a saddle point is the 
set of all regular points whose forward (backward) 
integral curve emerges from the saddle point. Sec¬ 
tion 2 contains a more precise definition. A non¬ 
degenerate saddle point has two stable and two un¬ 
stable separatrices. A Morse-Smale function is a 
Morse function whose stable and unstable separatri¬ 
ces are disjoint. In particular, the unstable separa¬ 
trices fiow into a sink, and separate the unstable re¬ 
gions of two sources. Similarly, the stable separatri¬ 
ces emerge from a source, and separate the stable re¬ 


gions of two sinks. The corresponding gradient vector 
field is called a Morse-Smale system (MS-system). 
The Morse-Smale complex (MS-complex for short) 
is a complex consisting of all singularities, separatri¬ 
ces and the open regions forming their complement, 
of the MS-system. In other words, a cell of the MS- 
complex is a maximal connected subset of the domain 
of h consisting of points whose integral curves have 
the same origin and destination. See also [10, 11, 22] 
and Section 2. The MS-complex describes the global 
structure of a Morse-Smale function. 

Existing algorithms for MS-complexes [10, 11] 
compute the complex of a piecewise linear function 
on a piecewise linear manifold, or, in other words, of 
a discrete gradient-like vector field. When h is an an¬ 
alytic function, we cannot use these algorithms with¬ 
out first creating a piecewise linear approximation h.. 
However, the MS-complex of h. is not guaranteed to 
be combinatorially equivalent to the MS-complex of 
the smooth vector field. The topological correctness 
depends on how close the approximation h is to h. 
Here “topological correctness” of the computed MS- 
complex M means that there is a homeomorphism f 
of the domain that induces a homeomorphism of each 
cell c G M to a cell c G M where M is the true MS- 
complex, and, moreover, this induced map c > c is 
an isomorphism of M and M. An isomorphism of two 
MS-complexes preserves the types of cells and their 
incidence relations. We can also require f to be an 
e-homeomorphism for some specified e > 0, i.e., the 
distance of a point and its f-image does not exceed 
£. As far as we know this problem has never been 
rigorously studied. Therefore, the main problem of 
this paper is to compute a piecewise-linear com¬ 
plex that is t-homeomorphic to the MS-complex 
of a smooth Morse-Smale function h. In short, we 
seek an exact computation in the sense of the Ex¬ 
act Geometric Computation (EGG) paradigm [17]. 
Note that it is unclear whether many fundamental 
problems from analysis are exactly computable in the 
EGG sense. In particular, the current state-of-the-art 
in EGG does not (yet) provide a good approach for 
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coping with degenerate situations, and, in fact, this 
paradigm needs to be extended to incorporate degen¬ 
eracies. Therefore, we have to assume that the gra¬ 
dients we start out with are Morse-Smale systems. 
However, generic gradients are Morse-Smale sys¬ 
tems [22], so the presence of degenerate singularities 
and of saddle-saddle connections is exceptional. Note 
that in restricted contexts, like the class of polyno¬ 
mial functions, absence of degenerate critical points 
(the first, and local, Morse-Smale condition) can be 
detected. However, even (most) polynomial gradient 
systems cannot be integrated explicitly, so absence 
of saddle-saddle connections (the second, and global, 
Morse-Smale condition) cannot be detected with cur¬ 
rent approaches. Detecting such connections even in 
a restricted context remains a challenging open prob¬ 
lem. 

1.2 Our contribution 

We present an algorithm for computing such a cer¬ 
tified approximation of the MS-complex of a given 
smooth Morse-Smale function on the plane, as illus¬ 
trated in Figures l(a)-(b). In particular, the algo¬ 
rithm produces: 

• (arbitrarily small) isolated certified boxes each 
containing a unique saddle, source or sink; 

• certified initial and terminal intervals (on the 
boundary of saddleboxes), each of which is guar¬ 
anteed to contain a unique point corresponding 
to a stable or unstable separatrix; 

• disjoint certified funnels (strips) around each 
separatrix, each of which contains exactly one 
separatrix and can be as close to the separatrix 
as desired. 

Note. The current version is an extensive elabora¬ 
tion of our previous paper [7] by incorporating all the 
theoretical results necessary to establish our method 
of certified Morse-Smale complex computation. The 


aim in [7] was more on providing an water-tight algo¬ 
rithm; however, the scope of showing all the theoret¬ 
ical details was limited. We complete that analysis 
part in the current extensive version. In Section 3, 
under certified critical-box computation, we provide 
the details of the relevant lemmas which were miss¬ 
ing in [7]. In Section 4, we establish rigorous theoret¬ 
ical foundations for refining the saddle-, source- and 
sink boxes that are used in computing the initial and 
terminating intervals of the stable and unstable sep- 
aratrices. All the theoretical results in this section 
are new additions to the current version. The final 
method section (Section 5) for the computation of 
disjoint certified funnels (strips) is now restructured 
into three subsections, each completes the relevant 
theoretical and algorithmic analysis. 

1.3 Overview 

Section 2 starts with a brief review of Morse-Smale 
systems, their singular points and their invariant 
manifolds. We also recall the basics of Interval Arith¬ 
metic, the computational context which provides us 
with the necessary certified methods. The construc¬ 
tion of the Morse-Smale complex of a gradient sys¬ 
tem Vh consists of two main steps: constructing dis¬ 
joint certified boxes for its singular points, and con¬ 
structing disjoint certified strips (funnels) enclosing 
its separatrices. Singular points of the gradient sys¬ 
tem are computed by solving the system of equations 
hx{x,y) = 0, h,y(x,y) = 0. This is a special instance 
of the more general problem of solving a generic sys¬ 
tem of two equations f(x,y) = 0, g(x,y) = 0. Generic 
means that the Jacobi matrix at any solution is non¬ 
singular, or, geometrically speaking, that the two 
curves f = 0 and g = 0 intersect transversally. In our 
context, this genericity condition reduces to the fact 
that at singular points of the gradient Vh. the Hessian 
matrix is non-singular. Section 3 presents a method 
to compute disjoint isolating boxes for all solutions 
of such generic systems of two equations in two un¬ 
knowns. This method yields disjoint isolating boxes 
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for the singular points of the gradient system. In Sec¬ 
tion 4 these boxes are refined further. Saddle-boxes 
are augmented with four disjoint intervals in their 
boundary, one for each intersection of the boundary 
with the stable and unstable separatrices of the en¬ 
closed saddle point. We also show that these inter¬ 
vals can be made arbitrarily small, which is crucial in 
the second stage of the algorithm. Sink- and source- 
boxes are refined by computing boxes—not necessar¬ 
ily axis-aligned—around the sink or source on the 
boundaries of which the gradient system is transver¬ 
sal (pointing into the sink-box and out of the source- 
box). This implies that all integral curves reaching 
(emerging from) such a refined sink-box (source-box) 
lie inside this box beyond (before) the point of inter¬ 
section. 

Section 5 describes the second stage of the algo¬ 
rithm, in which isolating strips (funnels) for the sta¬ 
ble and unstable separatrices are constructed. The 
boundary curves of funnels enclosing an unstable sep- 
aratrix are polylines with initial point on a saddle box 
and terminal point on a sink box. The gradient vec¬ 
tor field is transversally pointing inward at each point 
of these polylines. The initial points of the polylines 
are connected by the unstable interval through which 
the separatrix leaves its saddle box, and, hence, en¬ 
ters the funnel. The terminal points of these poly¬ 
lines lie on the boundary of the same sink-box. See 
also Figure 14. Given this direction of the gradient 
system on the boundary of the funnel, the unstable 
separatrix enters the sink-box and tends to the en¬ 
closed sink, which is its co-limit. Although the width 
of the funnel may grow exponentially in the distance 
from the saddle-box, this growth is controlled. We 
exploit the computable (although very conservative) 
upper bound on this growth rate to obtain funnels 
that isolate separatrices from each other, and, hence, 
form a good approximation of the Morse-Smale com¬ 
plex together with the source- and sink-boxes. These 
upper bounds are also used to prove that the algo¬ 
rithm, which may need several subdivision steps, ter¬ 
minates. 


We have implemented this algorithm, using Inter¬ 
val Arithmetic. Section 6 presents sample output of 
our algorithm. A contains guaranteed error bounds 
for the Euler method for solving ordinary differential 
equations, and B sketches a method for narrowing 
the separatrix intervals in the boundaries of the sad¬ 
dle boxes. 

1.4 Related Work 

Milnor [20] provides a basic set-up for Morse the¬ 
ory. The survey paper [5], focusing on geometrical- 
topological properties of real functions, gives an ex¬ 
cellent overview of recent works on MS-complexes. 
Originally, Morse theory was developed for smooth 
functions on smooth manifolds. Banchoff [4] intro¬ 
duced the equivalent definition of critical points on 
polyhedral surfaces. Many of the recent develop¬ 
ments on MS-complexes are based on this definition. 
A completely different discrete version of Morse the¬ 
ory is provided by Forman [12]. 

Different methods for computation. In the lit¬ 
erature there are two different method for comput¬ 
ing the Morse-Smale complexes: (a) boundary based 
approaches and (b) region based approaches. Bound¬ 
ary based methods compute boundaries of the cells of 
the MS-complex, i.e., the integral curves connecting 
a saddle to a source, or a saddle to a sink [28, 3, 11]. 
On the other hand, watershed algorithms for im¬ 
age segmentation are considered as region based ap¬ 
proaches [19]. Edelsbrunner et.al [11] computes the 
Morse-Smale complex of piecewise linear manifolds 
using a paradigm called Simulation of Differentia¬ 
bility. In higher dimensions they give an algorithm 
for computing Morse Smale complexes of piecewise 
linear 3-manifolds [10]. 

Morse-Smale complexes have also been applied in 
shape analysis and data simplification. Computing 
MS-complexes is strongly related to vector field visu¬ 
alization [14]. In a similar context, designing vector 
fields on surfaces has been studied for many graph¬ 
ics applications [30]. Cazals et.al. [6] applied discrete 
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Morse theory to molecular shape analysis. 

This paper contributes to the emerging area of 
Exact Numerical Algorithms for geometric problems 
[29]. Recent algorithms of this genre (e.g., [23, 18]) 
are numerical subdivision algorithms based on inter¬ 
val function evaluation and sign evaluation. 

2 Preliminaries 

In this section we briefly review the necessary math¬ 
ematical background on Morse functions, Morse- 
Smale systems, their singular points and their invari¬ 
ant manifolds. We also recall the basics of our com¬ 
putational model and Interval Arithmetic which are 
necessary for our certified computation algorithm. 

2.1 Mathematical Background 

Morse functions. A function h. : 2? C — > R 

is called a Morse function if all its critical points 
are non-degenerate. The Morse lemma [20] states 
that near a non-degenerate critical point a it is pos¬ 
sible to choose local co-ordinates x,y in which h. is 
expressed as h.{x,'y) = h(a) ± x^ Existence of 

these local co-ordinates implies that non-degenerate 
critical points are isolated. The number of minus 
signs is called the index ihla) of h. at a. Thus a 
two variable Morse function has three types of non¬ 
degenerate critical points: minima (index 0), saddles 
(index 1) and maxima (index 2). 

Integral curves. An integral curve x : I C R — > 2? 
passing through a point po on 2? is a unique maximal 
curve satisfying: x(t) = Vh,(x(t)), x(0) = po, for 

all t in the interval I. Integral curves corresponding 
to the gradient vector field of a smooth function h.: 
2? —> R have the following properties: 

1. Two integral curves are either disjoint or same. 

2. The integral curves cover all the points of V. 

3. The integral curves of the gradient vector field of 
h form a partition of V. 

4. The integral curve x(t) through a critical point po 


of h. is the constant curve x(t) = po. 

5. The integral curve x(t) through a regular point 
p of h. is injective, and if lim x(t) or lim x(t) 

t^OO t—^ —OO 

exists, it is a critical point of h. This implies inte¬ 
gral curves corresponding to gradient vector field are 
never closed curves. 

6. The function h. is strictly increasing along the in¬ 
tegral curve of a regular point of h. 

7. Integral curves are perpendicular to the regular 
level sets of h. 

Stable and unstable manifolds. Consider the in¬ 
tegral curve x(t) passing through a point p. If the 
limit lim x(t] exists, it is called the co-limit of p and 

t—^OO 

is denoted by cu(p). Similarly, lim x(t) is called 

t—> —OO 

the a-limit of p and is denoted by a(p) - again pro¬ 
vided this limit exists. The stable manifold of a sin¬ 
gular point p is the set W®(p) — [q £ V \ (X)(q) = p}. 
Similarly, the unstable manifold of a singular point 
p is the set W^(p) = {q g 2? j a(q) = p}. Here we 
note that both W®(p) and W’^(p) contain the singu¬ 
lar point p itself [15]. 

Now, the stable and unstable manifolds of a saddle 
point are 1-dimensional manifolds. A stable manifold 
of a saddle point consists of two integral curves con¬ 
verging to the saddle point. Each of these integral 
curves (not including the saddle point) are called the 
stable separatrices of the saddle point. Similarly, an 
unstable manifold of a saddle point consists of two 
integral curves diverging from the saddle point and 
each of these integral curves (not including the sad¬ 
dle point) are called the unstable separatrices of the 
saddle point. 

The Morse-Smale complex. A Morse function on 
2? is called a Morse-Smale (MS) function if its stable 
and unstable separatrices are disjoint. In particular, 
a Morse-Smale function on a two-dimensional domain 
has no integral curve connecting two saddle points, 
since in that case a stable separatrix of one of the 
saddle points would coincide with an unstable sep¬ 
aratrix of the other saddle point. The MS-complex 
associated with a MS-function h. on 2? is the subdi- 


5 



.2 Computational Model 



® Sink © Saddle O Source 

- stable Separatrix - Unstable Separatrix 

Figure 2: Morse-Smale complex 

vision of V formed by the connected components of 
the intersections W®(p) fl W^{q), where p, q range 
over all singular points of h. If 2? = K.^, then, ac¬ 
cording to the Quadrangle Lemma [11], each region 
of the MS-complex is a quadrangle with vertices of 
index 0,1,2,1, in this order around the region. 

Stability of equilibrium points. We note that a 
gradient vector field of a MS-function h. : 2? —> R 
can have three kinds of equilibria or singular points, 
namely, sinks (corresponding to maxima of h), sad¬ 
dles (saddles of h.) and sources (corresponding to 
minima of h). These singular points can be distin¬ 
guished based on the local behavior of the integral 
curves around those points. Locally, a sink has a 
neighborhood, which is a stable 2-manifold. Simi¬ 
larly, locally a source has a neighborhood, which is 
an unstable 2-manifold. Locally, a saddle has a sta¬ 
ble 1-manifold and an unstable 1-manifold crossing 
each other at the saddle point. A sink is called a sta¬ 
ble equilibrium point, where as a source or a saddle 
is called unstable equilibrium point. We note that, 
a source corresponding to a MS-function h. is a sink 
corresponding to the function —h. 


Dur computational model has two simple founda- 
ions: (1) BigFloat packages and (2) interval arith- 
netic (lA) [21]. Together, these are able to deliver ef- 
icient and guaranteed results in the implementation 
)f our algorithms. A BigFloat package is a software 
mplementation of exact ring (+, —, x) operations, 
livision by 2, and exact comparisons, over the set 
^ = {m2^ ; m, n G Z} of dyadic numbers. In prac- 
ice, we can use IEEE machine arithmetic as a filter 
or BigFloat computation to speed up many compu- 
ation. Range functions form a basic tool of lA: 
;iven any function F : —> R.”^, a range function 

□ F for F computes for each m- dimensional interval 
I (i.e., an m-box) an n-dimensional interval OF]!), 
such that F(I) C OF]!]. A range function is said to 
be convergent if the diameter of the output interval 
converges to 0 when the diameter of the input interval 
shrinks to 0. Convergent range functions exist for the 
basic operators and functions, so all range functions 
are assumed to be convergent. Moreover, we assume 
that the sign of functions can be evaluated exactly 
at dyadic numbers. All our boxes are dyadic boxes, 
meaning that their corners have dyadic coordinates. 

Interval implicit function theorem. To intro¬ 
duce a useful tool from lA, we recall some notation 
for interval matrices. An n x n interval matrix □ M 
is defined as 

□ M={M|Mq G DMij,!,) G{1 ...n}} 

Also note that we write 

0 ^ det □ M 

if there exists no matrix M G □ M such that det M = 
0 where det represents the determinant of the corre¬ 
sponding matrix. 

If I = Ix X ly is a 2D-interval (box) in R^, the 
interval Jacobian determinant □ gfl’yj is the 2x2 
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interval determinant given by 


11 9 (f) g) 

0(x,y) 


(I) 


ox 01) 

d| 9(,) nli,,, 

ox or) 


Proposition 2.1. (Interval Implicit Function 
Theorem, Snyder [26, 27]) Let F : R x R —> K. be 
a -map with components f and g. // I C R x R 
is a box for which 0 4 (I), then the system 

f(x,r)) =0, g{x,y) = 0 has at most one solution 
in I. 


Remark 2.1. In this paper, the domain D of h. is 
a finite union of axis-aligned dyadic boxes. Further¬ 
more, we (have to) assume that all stable and unsta¬ 
ble separatrices of the saddle points are transversal to 
the boundary. Computationally this means that any 
sufficiently close approximation of these separatrices 
is transversal to the boundary as well. □ 


3 Isolating boxes for singulari¬ 
ties of gradient fields 

As a first step towards the construction of the Morse- 
Smale complex of h. we construct disjoint isolating 
boxes for the singular points of Vh. To this end, we 
first show how to compute isolating boxes for the so¬ 
lutions of a generic system of two equations in two 
unknowns, which are confined to a bounded domain 
in the plane. This domain is a finite union of dyadic 
boxes. Applying this general method to the case in 
which the two equations are defined by the compo¬ 
nents of the gradient vector field Vh. we obtain iso¬ 
lating boxes for the singularities of this gradient field. 

3.1 Certified solutions of systems of 
equations 

We consider a system of equations 

f(x,y)=0, g{x,y)=0, (1) 


where f and g are -functions defined on a bounded 
axisparallel box D) C R2 with dyadic vertices. Fur¬ 
thermore, we assume that the system has only non¬ 
degenerate solutions, i.e., the Jacobian determinant 
is non-zero at a solution. In other words. 


g) 

a(x,y) 


^0 

(xcUo) 


( 2 ) 


for (xo,yo) satisfying (1). Geometrically, this means 
that the curves given by f(x,y) = 0 and g(x,y) = 0 
are regular near a point of intersection (xo,yo) and 
the intersection is transversal. We will denote these 
curves by Zf and Zg, respectively. Note that this 
condition is satisfied by Morse-Smale systems, since 
in that case f = h^ and g = hy, so the Jacobian 
determinant is precisely the Hessian determinant. 

Since the domain D of f and g is compact and we 
assume (2), the system (1) has finitely many solutions 
in D. Our goal is to construct a collection of axis- 
aligned boxes Bi,..., Bra and Bj,..., B(,^ such that 
(i) box Bt is concentric with and strictly contained in 
B(, (ii) the boxes B( are disjoint, (hi) each solution 
of (1) is contained in one of the boxes B|, and (iv) 
each box B( contains exactly one solution (contained 
inside the enclosed box Bi). The box pair (B|,B() 
is certified: Bt contains a solution, and B( provides 
positive clearance to other solutions. In fact, the se¬ 
quence of boxes will satisfy the following stronger 
conditions; See also Figure 3. 


1. The curves Zf and Zg each intersect the bound¬ 
ary of B^ transversally in two points. 

2. There are disjoint intervals If(f), iflg), I| (f) 
and (g) in the boundary of B| (in this order), 
where the first and the third interval each con¬ 
tain one point of intersection of Zf and 0B|, and 
the second and fourth interval each contain one 
point of intersection of Zg and 0B^. 


3. The (interval) Jacobian determinant of f and g 
does not vanish on B|, i.e.. 


0 ^ □ 


9(fig) 

9(x,y) 
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intervals in its boundary for Zf fl 0B^ and Zg fl 
The curves Zf and Zg intersect inside B^ iff these 
intervals are interleaved. 

3.2 Construction of certified box pairs 

We first subdivide the domain D into equal-sized 
boxes I (called grid-boxes), until all boxes satisfy cer¬ 
tain conditions to be introduced now. For a (square, 
axis-aligned) box I, let Np(I) be the box obtained 
by multiplying box I from its center by a factor of 
1 + p, where j < P < 1. We also denote N i (I) by 
N (I) ; this is the box formed by the union of I and 
its eight neighbor grid-boxes. We shall call N (I) the 
surrounding box of I. The algorithm subdivides D 
until all grid-boxes I satisfy Co (I) V ( Ci (I) A C 2 (I)), 
where the clauses Ci(I), i = 0,1,2, are the following 
predicates: 

Co(I): 0 ^ Df)!) VO ^ Dg)!) 

Ci(I): 0^n|J^(N(I)) 

9(vy) 

C 2 (I): C 2 (I,f)AC 2 (I,g), 

where 

C 2 (I,f) = (□^(N(I)), □^(N(I))) > cos 

If Co(I) holds, box I does not contain a solution, 
so it is discarded. The second predicate guarantees 
that N(I) contains at most one solution. This is a 
consequence of Interval Implicit Function Theorem 
[25, 24]. Condition C 2 (I) is a small angle variation 


condition, guaranteeing that the variation of the unit 
normals of the curves Zf fl N (I) and Zg fl N (I) do not 
vary too much, so these curves are regular, and even 
‘nearly linear’ (the unit normal of Zf is the normal¬ 
ized gradient of f). Here (•,•) denotes the interval 
version of the standard inner product on K^. The 
lower bound for the angle variation of Vf is gener¬ 
ated by the proof of Lemma 3.2. 

Remark 3.1. Condition Ci (I) implies that there is 
a computable positive lower bound a(I) on the an¬ 
gle between Vf(p) and Vg{q) where p, q range over 
the surrounding box of I. More precisely, to com¬ 
pute a(I), we first compute a lower bound L on the 
quantity ||vf||.^||vg|| • ^his L may be obtained 

by an interval evaluation of this quantity at N(I); 
note that L > 0 iff condition Ci (I) holds. We define 
a(I) as arcsin(L). □ 

Our algorithm will construct disjoint certified 
boxes surrounding a box I. As observed earlier, the 
surrounding boxes N(I) and N(J) of disjoint boxes 
I and J may intersect. Since our algorithm will con¬ 
struct disjoint certified boxes surrounding a box I, its 
surrounding box should be smaller than the box N (I). 
To achieve this, note Np(I) be the box obtained by 
multiplying box I from its center by a factor of 1 + p, 
where y < p < 1. In particular, N (I) = N 1 (I). If 
N (JJ) nl = 0, then / 2 (I] and N j/ 2 (J) have disjoint 
interiors. This is a key observation with regard to 
the correctness of our algorithm. 

Lemma 3.2. Let I be a box such that conditions 
“■Coll), Ci(I), and C 2 (I) hold. Let d be the length 
of its edges, and let j < P< 1. 

^1. If Zf intersects I, it intersects the boundary 

^^/Np(I) transversally at exactly two points. At a 
point of intersection of Zf and an edge e o/0Np(I) 
the angle between Zf and e is at least 
2. If I contains a point of intersection of Zf and 
Zg, then the points of intersection of Zf and 
0Np(I) are at distance at least 2pdtan2«-(I) from 
the points of intersection of Zg and 0Np(I). On 
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9Np(I), the points of intersection with Zi and with 
Zg are alternating. 
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Figure 4: Lower bound on angle of intersection of Zf Figure 5: A lower bound for the distance between 
and the boundary of the surrounding box Np(I) in points of Zf and Zg on the boundary of Np(I). 
case Zf intersects I. 


Proof. 1. Assume that Zf intersects a vertical edge of 
9Np(I) at Let p be a point on Zf fll. See Figure 4. 
Then there is a point s on the curve segment at 
which the gradient of f is perpendicular to this line 
segment p^. Let |3 be the (smallest) angle between 
Vf(s) and the horizontal direction. Referring to the 
rightmost picture in Figure 4, we see that this angle 
is not less than Zpqr = jTt—Zqpr. Since || q— pi| = 
di/l + 2p + 2p^, it follows that 



P 

arccos — , 

Vl +2p + 2p2 


> 


7t 
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where the last inequality holds since j < P < 1. 
Since condition C 2 (N(I),f) holds, the angle between 
the gradients of f at s and q' is less than ^. There¬ 
fore, the angle between Zf and the vertical edge of 
9Np(I) at ■q is at least 

2. Let p G I be the point of intersection of Zf and 
Zg. Suppose Zf and Zg intersect the vertical edge e 
of 9Np(I) in q and r, respectively. Then there is a 
point s on Zf between p and q at which the gradient 
of f is perpendicular to pq, and there is a point on 
Zg between p and r at which the gradient of g is 
perpendicular to pr. Let a(I) be the lower bound on 
the angle between Vf(p) and Vg(q) where p, q range 


over the box I. See Figure 5, where a > a{I). For 
fixed a, the distance between q and r is minimal if 
the projection p of p on the edge e is the midpoint of 
qr, in which case |jr — q j] = 2||p — p|| tan^a. Since 
IIP “ P II > pd. and a > a(I), the distance between q 
and r is at least 2pdtan □ 

The following result gives an estimate on the posi¬ 
tion of the points at which Zf intersects the boundary 
of the surrounding box of I in case Zf intersects an 
edge of 91 in at least two points. For an edge e of the 
inner box I let I and r be the points of intersection 
of the line through e and the edges of the surround¬ 
ing box Np{I), perpendicular to e. See also Figure 6. 
The dyadic intervals on the boundary of this sur¬ 
rounding box with length at most :^(1 + p)d, cen¬ 
tered at I and r, respectively, are denoted by Lp(e) 
and Rp(e), where d is the length of the edges of I, 
and j < P < 1. Here, we note that in the foilwing 
lemma 3.3 intervals Lp(e) and Rp(e) can be made 
corresponding dyadic intervals by replacing its real 
endpoints with suitable conservative dyadic numbers 
satisfying the conditions of the lemma. 

Lemma 3.3. Let I be a box such that-'Co(I), Ci (I) 
and Cl (I) hold, and let e be one of its edges. Let 
I < P< 1. 
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1. If Zf intersects an edge e of the boundary of 
I in at least two points, then it transversally 
intersects 9Np(I) in exactly two points, one in 
each of the dyadic intervals Lp(e) and Rp(e). 

2. If Zf intersects 9Np(I) in the dyadic inter¬ 
vals Lp(e) and Rp(e), then these intersections 
are transversal, and Zf intersects 9Np(I) at 
exactly two points, one in each of these in¬ 
tervals. 



on Zf at which Vf is perpendicular to qs. Since the 
angle of qs and the vertical direction is at most 

^ (1 + P)/8 ^ 1 ^ ^ 

arctan—;-= arctan- < —, 

1 + p 8 25 ’ 

it follows from the small normal variation condition 
C 2 (I) that the gradient of f at any point of Np(I) 
makes an angle of at most ^ with the 

vertical direction. This rules out multiple intersec¬ 
tions with the vertical edges of 9Np(I). It also im¬ 
plies that Zf lies above the polyline qms, where m 
is the intersection of the line through q with slope 
pj^an-^ and the line through s with slope — tan-^. 
Therefore, all points of Zf lie at distance at most 
^(1 + p)d + ^(1 + p)d tan-p^ < pd from the line 
through e, so Zf does not intersect the edges of 
9Np(I) parallel to e. □ 


Figure 6: Intervals containing points Zf fl 9Np(I). 


Proof. 1. There is a point on Zf between the points 
of intersection with e at which the gradient of f is 
perpendicular to e. See Figure 6. The small nor¬ 
mal variation condition C 2 (N{I),f) implies that Zf 
does not intersect any of the two edges of 9Np(I) 
parallel to e, and that it intersects each of the edges 
of 9Np(I) perpendicular to e transversally in exactly 
one point. Let q be the point of intersection with 
the edge containing r. Then there is a point on Zf at 
which the gradient of f is perpendicular to the line 
segment pq. Since the angle between the gradients 
of f at two points of Np{I) does not differ by more 
than 3^71, we have Zqpt < ^7i. Therefore, 


3.3 Towards an algorithm 

After the first subdivision step, we have constructed 
a finite set B of boxes, all of the same size, such that 
Co (I) A Cl (N (I)) A C 2 (N(I)) holds for each box I. 
For each grid-box I, the algorithm calls one of the 
following: 

• Discard(I), if it decides that I does not contain 
a solution. It marks box I as processed. 

• ReportSolution(I). It returns the certified 

pair (Ni /2 (I)> and marks all boxes con¬ 

tained in N (I) as processed. 

In the latter case a solution is found inside N i /2 (I), 
but, as will become clear later, it may not be con¬ 
tained in the smaller box I. In view of Ci (I) none 
of the grid-boxes in N (I) contain a solution different 


q—r|| = II p—TII tanZqpr < (l+p)dtan^ < -^(l+plflom the one reported, so they are marked as being 

processed. 

Decisions are based on evaluation of the signs of f 
and g at the vertices of the grid-boxes (or at certain 
dyadic points on edges of grid-boxes). An edge of a 
box is called bichromatic {monochromatic) for f if 
the signs of the value of f at its vertices are opposite 
(equal). 


In other words, Zf intersects Rp. Similarly, Zf in¬ 
tersects Lp. The small normal variation condition 
C 2 (N (I) ) implies that there are no other intersections 
with the edges of 9Np(I). 

2. Let the points of intersection of Zf and Rp(e) and 
Lp(e) be q and s, respectively. Then there is a point 
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Algorithm, case 1: I has a bichromatic edge 
for f and a bichromatic edge for g Then Zf 
and Zg intersect I, and, according to Lemma 3.2, 
part 1 , both curves intersect the boundary of 72(1) 
transversally in exactly two points. For each of the 
two points in 9 Ni/ 2 (I) the algorithm computes an 
isolating interval—called an f-interval—on 0N i 72(1) 
of length 2 dtan 2 a(I). The two g-intervals are com¬ 
puted similarly. If the f- and g-intervals are not 
interleaving, there is no solution of ( 1 ) in box I— 
even though there may be a solution in N ] 72 ( 1 )—and 
Discard(I) is called. This follows from Lemma 3.2, 
part 2. If the intervals are interleaving, then there is a 
point of intersection inside N ] /2 (I) , so the algorithm 
calls ReportSolution(I). 

Algorithm, case 2: I contains no bichromatic 
edge for f (g), and at least one bichromatic 
edge for g (f, respectively) We only consider 
the case in which all edges of I are monochromatic 
for f. Then the algorithm also evaluates the sign of f 
at the vertices of the box N 1 /2 (I) • If N 172 (I) has no 
disjoint bichromatic edges (as in the fourth and fifth 
configuration of Figure 7), the isocurve Zf does not 
intersect I, so the algorithm calls Discard(I). To 
deal with the remaining case, in which N 1 /2 (I) has 
two disjoint bichromatic edges (as in the sixth con¬ 
figuration in Figure 7) we need to evaluate the sign of 
f at certain dyadic points of these bichromatic edges, 
followed from Lemma 3.3. By evaluating the signs 
of f at the (dyadic) endpoints of the interval Lp(e) 
and Rp(e) the algorithm decides whether they con¬ 
tain a point of intersection with Zf. If at least one of 
these intervals is disjoint from Zf, then Discard(I) 
is called. Otherwise, the algorithm computes isolat¬ 
ing f- and g-intervals of length j d tan j a(I). As in 
case 1, the algorithm calls ReportSolution(I) if 
these intervals are interleaving, and Discard(I) oth¬ 
erwise. 

Algorithm, case 3: all edges of I are monochro¬ 
matic for both f and g Again, let e be the 
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Figure 7: Sign patterns of the box Ni/ 2 (I) enclosing 
the grid-box I with monochromatic edges for f. The 
three top configurations are ruled out by the small 
normal variation condition C 2 (I). The fourth, fifth 
and sixth configuration are all possible, but only in 
the sixth situation Zf may intersect the inner box. 


(unique) edge of I closest to the edge of N ] /2 (I) which 
is monochromatic for f, at whose vertices the sign of 
f is the opposite of the sign of f at the vertices of I. 
Edge e' of I is defined similarly for g. 

Case 3.1: e = e'. In this case Zf or Zg does not 
intersect I. Indeed, if Zf intersects I, it intersects e 
in at least two points, so there is a point p G Zf at 
which Vf(p) is perpendicular to e. Condition Ci (I) 
guarantees that Vg is nowhere parallel to Vf(p), so 
Zg does not intersect e, and, hence, does not inter¬ 
sect I. Therefore, Discard(I) is called. 

Case 3.2: e ^ e'. If Zf does not intersect 
Lp(e) or Rp(e), or if Zg does not intersect Lp(e') 
or Rp(e'), then, as in case 2, the algorithm calls 
Discard(I). Otherwise, Lp(e) or Rp(e) are isolat¬ 
ing f-intervals which are disjoint from the isolating 
g-intervals Lp(e') or Rp(e'). If e and e' are per¬ 
pendicular, then these f- and g-intervals are inter¬ 
leaving, and, hence, ReportSolution(I) is called. 
Otherwise, there is no solution in I, so Discard(I) 
is called. 
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Refinement: disjoint surrounding boxes We 
would like distinct isolating boxes I, JJ to have disjoint 
surrounding boxes N (I), N (JJ). There is a simple way 
to ensure this: we just use the predicate Ci(N(I)) 
to instead of Ci (I) in the above subdivision process. 
Then, if the interior of N (I) fl N (J) is non-empty, we 
can discard any one of I or JJ. 


4 Isolating boxes for sinks, 
sources and saddles 

In a first step, described in Section 3, we have con¬ 
structed certified disjoint isolating boxes ...,B,^ 
the singular points of Vh. in the domain JD of h.. Let 
D* be the closure of D \ (JB-, U • • • U 1B„^). Obviously, 
D* is a compact subset of . 

In a second step towards the construction of the 
MS-complex, we refine the saddle-, sink- and source- 
boxes. In Section 4.1 we show how to augment each 
saddlebox by computing four arbitrarily small dis¬ 
joint intervals in its boundary, one for each intersec¬ 
tion of a stable or unstable separatrix with the box 
boundary. Subsequently, in Section 4.2, we show how 
to construct for each source or sink of Vh (minimum 
or maximum of h) a box on the boundary of which 
the gradient field is pointing outward or inward, re¬ 
spectively. These boxes are contained in the source- 
and sinkboxes constructed in the previous section, 
but are not necessarily axes-aligned. 


4.1 Refining saddle boxes: Isolating 
separatrix intervals 

To compute disjoint certified separatrix intervals we 
consider wedge shaped regions with apex at the sad¬ 
dle point, enclosing the unstable and stable manifolds 
of the saddle point. Even though the saddle point is 
not known exactly, we will show how to determine 
certified intervals for the intersection of these wedges 
and the boundary of a saddle box. 


First we determine the eigenvalues and eigenvec¬ 
tors of the Hessian of h at a point (xo,yo) in the 
interior of the saddlebox I —its center point, say— 
and consider these as good approximations to the 
eigenvalues and eigenvectors of the Hessian, i.e., the 
linear part of Vh., at the saddle point. Let H be the 
Hessian, i.e.. 


H = 




(3) 


and let be the Hessian evaluated at (xo,yo)- The 
eigenvalues Au and As of are given by 


Au = I (h°^ -h h°y -h iy(h0,,-h0y)2 -h4{h0y)2) 
As = I (H°,, -h h°y - ^(h0,,-h0y)2 -h4{h0y)2), 


and the corresponding eigenvectors are 



The singular point is a saddle, so we have Aj < 0 < 
Au. Since is a symmetric matrix, its eigenvectors 
are orthogonal. More precisely. 


V® = 



V2(V-). 


Here denotes counterclockwise rotation over an 
angle a. Therefore, 

IIV'^II = ||V"^!|, and det(V^,V'^) = ||V"^!|^ (5) 


The stable and unstable eigenvectors V® and are 
good approximations of the tangent vectors of the 
stable and unstable manifolds of the saddle point. 
These invariant manifolds are contained in wedge- 
shaped regions, which are defined as follows. 

Definition 4.1. Let the (orthogonal) vectors 
and V® he the stable and unstable eigenvectors of 
the Hessian of h at the center of a saddle box 
I, and let |3 G (0,^). The unstable wedge Cp is 
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the set of points in the surrounding box N (I) at 
which the (unsigned) angle between Vh. and 
is at most |3. See Figure 8. Similarly, the stable 
wedge Cp is the set of points in the surrounding 
box N(I) at which the (unsigned) angle between 
Vh. and V® is at most |3. 

'“P 


N(I) 

Figure 8 : The unstable wedge Cp enclosing the un¬ 
stable separatrix, and the stable wedge Cp enclosing 
the stable separatrix. The gradient vector field Vh, 
represented by solid arrows, is transversally pointing 
inward along the boundary of the unstable wedge, 
and outward along the boundary of the stable wedge. 
At points of the unstable wedge boundary U Fhp 
the vector field Xp is parallel to or — V’^, so Vh 
makes an angle —(3 with or — there. 

The saddle point belongs to both the stable and the 
unstable wedge. Since and V® are orthogonal, 
and 0 < (3 < ^, this is the only common point of the 
stable and unstable wedge. 

Conditions We now introduce additional condi¬ 
tions, which guarantee that each of the wedge- 
boundaries consist of two regular curves, cf 
lemma 4.2. In fact, these conditions guarantee that 
Cp and Cp are really wedge-shaped. Fix a > 1, 
and let 5 > 0 be an arbitrarily small constant (to 


be specified later). At the point (xo,r)o) we have 
HV^ = AuV^, HV® = AsV% so N(I) can be taken 
small enough to guarantee that the following condi¬ 
tion is satisfied at all points of N (I): 

Condition I(a,I). At every point of the box N(I) 
the following inequalities hold: 

-Au ■ IIV'" II < II II < aA^, • II||, 

a 

-|As|-||V^||<||HV*||<a|As|-||V*||. 

a 

At the point (xo,yo) we also have = 

AullV^^lP, (HV^V'^) = AsllV^lP, and (HV%V^) = 

0. Therefore, for any 6 > 0, the box N(I) can be 
taken small enough such that the following condition 
holds: 

Condition 11(5,1). At every point of the box N(I) 


the following inequalities hold: 


(HV^,V^) > ^A^,||V ^||2 

( 6 ) 

(HV^V^) > IIA 3 I ||V -||2 

(7) 

|(HV%V^)| < 6||W"||2 

( 8 ) 


Since H is symmetric, ( 8 ) also implies V®)| < 

5|| V^lp. 


On the boundary of Cp the gradient field makes a 
(signed) angle ±(3 or tt ± (3 with or, in other 
words, X±p is (anti)parallel to Again, Xp is the 
vector field obtained by rotating Vh. over an angle 
|3. So let r^p be the curve along which the vector 
field X ±(3 is (anti)parallel to the unstable eigenvector 
V’^. Then the boundary of the unstable wedge is the 
union of the two curves Fp and FF^p. The curve F|^ is 
defined by the equation 

det(V"",Xp(x,'y)) = 0. (9) 

Obviously, the saddle point lies on Fp. The function 
rjj’^p is defined similarly. 

Similarly, the boundary of the stable wedge is the 
union of curves F|.p, along which the vector field X ±(5 
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is (anti)parallel to V®. The curves r|.p are defined by so 
the equation 

t|)±l 3 (x,y) := det(V*,X±| 3 (x,'y)) =0. 


The following technical result provides computable 
upper bounds for the angle variation of the normals 
of the boundary curves of the stable and unstable 
wedges. 

Lemma 4.2. Let coi G (0, (to be specified 
later), let a > 1, and let I be such that Condi¬ 
tion l(a,I) holds. Let 0 < |3 < j and 6 > 0 such 
that 


1 = 

bJ 

Therefore, 


cos |3 sin|3\ /—VJ 
— sin |3 cos |3 / \ 


= (cos|3)V*+(sin|3) V^. 




^hxx hxy \ fA 


) = COS |3 (HV'^j+sin |3 (HV'"). 
'13/ 


. sin 0)1 . I As I I Auk 

sm|3 < -—- 7 = mm( — , — j, 


(15) 

Condition I (a, I) implies that 0 and HV® 

0, and {HV^, HV®} are independent vectors, at all 
points of N(I), so the gradient of ijjp is nonzero at 
every point of N (I), so is a regular curve. 
Expression (15) for ViJjp implies that 


6 < 

5 < 


4a^\/2 

sin cui 


8 a 
tan |3 


min(|As|, |Au|), 


(10) II 11^ = cos^ |3 II HV® ||^+2sin |3 cos |3 (HV% HV"")+sin^ |3 || HV^ 
Usin 

( 11 ) 


min(|As|, |Au|). 


( 12 ) 


If Condition 1I(S,I) also holds, then at any point 
o/N(I) 


Using the Cauchy-Schwarz inequality 
|(HV%HV^)| < IIHVMI • IIHV^'II and the fact 
that |3 > 0 we get 

II Vrl)p 11^ > cos^ |3 II HV® - 2 sin |3 cos |3 || HV® || • || HV 
+ sin^ |3 II HV^ 11^ 

= cos^ (3 (II HV® II - IIHV^ II tan |3)^ (16) 


71 


7t 


7t 


-- 0)1 < angle(Vr|)e,V'^) < < angle(Vij)p, V^) < -j.+o)!.. sinuJi , A 


P’ ' ’ Since'sin |3< , 

(13) 402^2'A 


-— and 0 < |3 < ?, it follows 


and 


from Condition l(a, I) that 

7t 


2 0)1 < angle)Vr|)p,V®) < ^ < angle(Vrl)ip, V®] < || tan|3 < oAu II|| < — || V''^||. 

(14) 

j. j. Using Condition I again we get 

In particular, the angle variation of any of the 

gradients Vi|)^p and Vi|)^p over N (I) is less than || [[ || || tan |3 > ||yu|| || || V^H 


2o) 


Proof. We only show that the angle variation of Viji p 
over N(I) is less than 2cui. Since 


2 a 

In view of (16) we get, using cos |3 > 


2a 


Xp =h,. 


(cos |3 
(sin |3 


ivii);^ll> 


|As 


iv^ll. 


+ h-ii 


^ 2 aV 2 

Expression (15) for Vijip also implies that 
the function i|)p satisfies i()p = Ap hx + Bp hy , where = cos |3 (HV®, V^) + sin |3 (HV^, V^) 


(17) 


— sin |3\ 
cos |3 / ’ 


Ap = det(V-, 


cos |3 


sin |3 


and Bp = det(V^, 


— sin |3^ 
cos |3 


cosP(HV-,V-){1^^^7^+tanP). 
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Condition II and (12) imply 




26 1 

< — < j tan |3. 
Au 


Since |3 > 0, this implies (Vx|)p, V’^) > 0 on N(I). 

According to Condition II we have |(HV®,V’^)| < 
6 II Ip, whereas the Cauchy-Schwarz inequality im¬ 
plies |(HV"^,V^)| < • IIV^II < aAullV^lp. 

Therefore, 


Corollary 4.3. Let a and cui be constants such 
that a > 1, and cui = | arctan^. Let |3 G (0, cui) 
and 6 > 0 such that (10), (11) and (12) hold. 

//1 is a saddle box with concentric surrounding 
box N(I) satisfying Condition I (a, I) and Condi¬ 
tion 11(6,1), then 

1. The saddle point is the only common point of 
the stable wedge and the unstable wedge 

ps 

'“P- 


0 < V^) < (6 cos |3 + aAu sin |3) || H^. 

Together with (17) this implies 

2aV2 

lv^> < (6cos 13 + sin W 

Given the upper bounds (11) for 6 and (10) for sin |3, 
we get 

2aV2 o n 4a ^ , laVl , . 

——— 6 cos p < ■— 6 < j sm 0)1 and aA^ sm 

|Asl |As| |As| 

Therefore, 

0< ( ||y^u|| »p^) < =cos(f-a)i) 

At all points of N (I) we then have 

^ - 0)1 < angle(Vxj)p,V^) < 

Since is constant, the angle variation of Vt|)p over 
N(I) does not exceed 2o)i. □ 

The main result of this subsection states that, 
under suitable conditions, the intersection of the 
boundary of a saddle box and the stable and un¬ 
stable wedges can be computed. Moreover, at all 
points of these intersections the gradient vector field 
is transversal to the boundary of the saddle box, and, 
even stronger, at these points there is a computable 
positive lower bound for the angle of the gradient 
vector field and the boundary of the saddle box. 


2. The gradient vector field Vh is transversal at 
points on the boundary of these wedges, dif¬ 
ferent from the saddle point: on the boundary 
of the unstable wedge it points inward, except 
at the saddle point, and on the boundary of 
the stable wedge it points outward, except at 
the saddle point. 

3. The unstable wedge contains the unstable 
separatrices of the saddle point, and the stable 

^ ]We^^ Cp contains the stable separatrices. 

4 . The unstable wedge intersects the boundary of 
N(I) in two intervals, called the unstable in¬ 
tervals. Similarly, the stable wedge intersects 
the boundary of N (I) in two intervals, called 
the stable intervals. These four intervals are 
disjoint, and the unstable and stable intervals 
occur altematingly on the boundary of N (I). 
At each point of a stable or unstable inter¬ 
val the (unsigned) angle between Vh and the 
boundary edge containing this point is at least 
0 ) 1 . Moreover, there are computable isolating 
intervals for each stable and unstable inter¬ 
val. 

Proof. 1. If p G Cp n Cp, then Vh(p) makes an 
angle |3 G (0,^) with both and V®. There¬ 
fore, Vh(p) = 0, since these vectors are orthogonal. 
Hence, p is a singular point of Vh inside N (I), which 
is the saddle point. 

2. Recall from Lemma 4.2 that (Vx|)p, V^) is positive 
on N (I). Let p be the saddle point of Vh in I. Since 
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Xp is parallel to on one component of Tp \{p}, and 
parallel to —on the other component, it follows 
that Xp is pointing into the unstable wedge along 
both components. See again Figure 8. Since Vh. is 
obtained by rotating Xp over —13, and X_p over |3, 
also VH is pointing into the unstable wedge. Simi¬ 
larly, (Vx[>’^p, is negative along both components 
of r^p \ {p}, so VH is also pointing into the unstable 
wedge along each of these components. 

A similar argument shows that VH is pointing out¬ 
ward along each of the boundary components of the 
stable wedge Cp, except at the saddle point. 

3. The second part of the lemma implies that the un¬ 
stable wedge is forward invariant under the flow of 
the gradient vector fleld VH. In particular, it con¬ 
tains the unstable separatrices of the saddle point. 
Similarly, the stable wedge Cp is backward invari¬ 
ant, so it contains the stable separatrices of the sad¬ 
dle point. 

4. Suppose Pp intersects an edge e of the surrounding 



Figure 9: Lower bound on angle of intersection of Pp 
and the boundary of the surrounding box N (I). 


box at a point 'q, see Figure 9. We first show that the 
angle of VH('q) and e is bounded away from zero. To 
see this, observe that there is a point s G Pp at which 
Vi|)p (s) is perpendicular to p'q. Therefore, the angle 
between Vt|)p{s) and the normal of e is at least cuo, 
where cuo = arctan j — 3a)i. 

The angle between and ViJ)p lies in the interval 
— o)], ^), cf Lemma 4.2, so the angle between 


and e is at least cuq — cui = 2a)i. 

At any point of Cp the angle between VH and is 
at most |3 —by the definition of Cp—so at any point 
of G Cp n e the angle between VH and e is at least 
2a)i — (3 > 0)1. 

To And isolating intervals for the intersection of the 
stable and unstable wedges with the boundary of the 
surrounding box N (I) , we compute isolating intervals 
for the intersection of each of the four curves 4’±p = 
0,ij)^p = 0 with this boudary. The normal to each 
of the curves 4’±p = 0 makes an angle of at least j — 
4a) 1 with each of the curves 4’±p = 0. This follows 
from (13) and (14), and the fact that and V® are 
perpendicular. Since cui = -j arctan^ < ^, so the 
angle between each of the curves i])^p = 0 and each 
of the curves 4’±p = 0 is at least which 

is bounded away from zero. Therefore, the method 
of Section 3 provides such isolating intervals. □ 

As will become clear in the certified construction 
of the MS-complex, we need to be able to provide cer¬ 
tified separatrix intervals of arbitrarily small width, 
without refining the saddle box: 

Lemma 4.4. Let S be a separatrix box satisfying 
the conditions of Corollary 4-3. Then the isolat¬ 
ing separatrix intervals in the boundary of N(I) 
can be made arbitrarily small. 

The proof of this result is rather technical. For a 
sketch we refer to B. 

4.2 Refining boxes for mcixinia and 
minima 

To construct the MS-complex, the algorithm needs 
to determine when an unstable (stable) separatrix 
will have a given maximum (minimum) of H as its 
a)-limit (a- limit). For each maximum (minimum) 
the algorithm determines a certified box such that 
the gradient vector fleld points inward (outward) on 
the boundary of the box. Unfortunately, we cannot 
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always choose an axis-aligned box, as will become 
clear from the following example. 

Let h.(x,y) = —5x^ — 4xy — then the origin is 
a sink of the gradient vector field 

This vector field is horizontal along the line y = “2x, 
which intersects the horizontal edges of every axis 
aligned box centered at the sink (0,0). In other 
words, the vector field is not transversal on the 
boundary of any such box. 

However, there is a box aligned with (approxima¬ 
tions) of a pair of eigenvectors of the linear part of the 
gradient vector field at (or, near) its singular point, 
for which the vector field is transversal to the bound¬ 
ary. To see this, we refine the sink-box to obtain three 
concentric axis aligned boxes Ji C Ji C JJ3, such that 

(i) the edge length of Ji, i = 2,3, is three times the 
edge length of Ji_i, and 

(ii) the sink is contained in the inner box Ji. See also 
Figure 10. Moreover, let Vi and V 2 be the (orthog- 



Figure 10: Construction of a sinkbox. 

onal) eigenvectors of the Hessian matrix at the 
center of the boxes. These eigenvectors, correspond¬ 
ing to the eigenvectors Ai and Ai, are computed as 
in Section 4.1, cf (4). 

We require that 


(iii) the gradients of the two functions t))] and t)) 2 , 
defined by 

^l^i(x,y) = (Vh(x,y), Vi), i = 1,2, 

have small angle variation over the outer box J3. This 
condition is made precise in Lemma 4.5 below. Note 
that 

Vrl)i(x,y) = H(x,y)Vi, (18) 

where H(x,y) is again the Hessian matrix of h. at 
(x,y). Since this matrix is non-singular, we can find 
a triple of boxes Ji C J2 C J3, satisfying conditions 
(i), (ii) and (iii), such that H(x,y) is nearly constant 
over the outer box (again, this is made precise in 
Lemma 4.5). In particular, Vi(>i is nearly parallel to 
Vi, since H°Vi = AiVi. Now construct boxes Ii and 
I 2 , which are the smallest boxes enclosing Ji and J2, 
respectively, with edges parallel to V^ or V^. 

Lemma 4.5. Suppose on the outer box J3 the fol¬ 
lowing conditions hold: 

1- l|HVi|!>i|Ai!.||Vi||,/or i = 1,2; 

2. |(HVi,V 2)| = |(HV2,Vi)| < 2||Vi Iparctani. 

Then the gradient vector field is transversal to the 
boundary of I 2 ■ 

Proof. The second condition limits the variation of 
the angle of Vi(>i = HVi and the basis vectors V] and 
V 2 over J3. Using this bound, we use the same ar¬ 
guments as in Section 3, applied to the pair of boxes 
Ii , I2, to show that the curve \l)i(x,y) =0 does not 
intersect the edges of I 2 perpendicular to Vi. There¬ 
fore, i(>i = (Vh, Vi) is nowhere zero on these edges, 
so Vh. is nowhere tangent to these edges. Therefore, 
I 2 is the desired sink box, on the boundary of which 
Vh is pointing inward. In other words, if an unsta¬ 
ble separatrix intersects the boundary of this box, the 
part of the separatrix beyond this point of intersec¬ 
tion lies inside the sink-box. Certified source-boxes 
are constructed similarly. □ 
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5 Isolating funnels around sep- 
aratrices 

If the forward orbits of the endpoints of an unsta¬ 
ble separatrix interval have the same sink of Vh. as 
co-limit, these forward orbits bound a region around 
the unstable separatrix leaving the saddle box via 
this unstable segment. This region is called a fun¬ 
nel for the separatrix (this terminology is borrowed 
from [16]). 

In this section we provide the details of the con¬ 
struction of a certified funnel for each separatrix. 
First we show, in Section 5.1 how to construct two 
polylines per separatrix interval those are candidates 
for the funnel around the corresponding separatrix. 
Then, in Section 5.2, we introduce the notion of 
width of a funnel, and derive upper bounds for its 
growth. These upper bounds are the ingredients for a 
certified algorithm computing these funnels. The al¬ 
gorithm computes the Morse-Smale complex by pro¬ 
viding disjoint certified funnels for each stable and 
unstable separatrix. The proof of correctness and 
termination is presented in Section 5.3. 

5.1 Construction of fences around a 
separatrix 

Let be the vector field obtained by rotating the 
vector field X = Vh. over an angle i.e., 

/cosfi — sinfi\ /h,x\ 

® ^ ^vsin^ cos^ ) IvHyj ’ 

We compute an isolating funnel for the forward or¬ 
bit of Vh. through a point p by enclosing it between 
(approximations of) the forward orbits of Xa and X_^ 
through p. See Figure 11. 

Small angle variation We first determine some 
bounds on the angle variation of the gradient vector 
field Vh over D*. We subdivide the region D* into 
square boxes over which the angle variation of Vh is 
at most fi, where ^ is to be determined later. Let 



Figure 11: Orbits of the rotated vector fields X^ and 
X_^ through a point p enclose the forward orbit of 
Vh through p. On the right polygonal lines approx¬ 
imating these orbits. 


w be the edge length of the boxes. If X = (f, g) is 
a vector field on then the angle variation over a 
regular curve P is given by [2, Section 36.7]: 

gdf-fdg 

Jr f^ + g^ • 

If X = Vh, this angle variation is equal to 

(hxhx^ - h^hxx) dx + (hxh^^ - h^h^^) dg 
. r + 

Let Co and Ci be constants such that 


max 


’chxy hyhx 


hi-hh2 


< Co and max 


hxh^^ - h^l 


h2-hh2 

(19) 

Then the angle variation over a curve P is less than 


J (Co dx + Cl dy) < (Co + Ci) length(P). 

This inequality provides an upper bound for the max¬ 
imal angle variation over a square box: 

Lemma 5.1. Let Co and Ci satisfy (19). Then 
the total angle variation over a square box in D* 
with edge length w does not exceed (Co + Ci )w-\/2. 

The grid boxes have edge length w such that the 
angle variation of Vh over any box in D* is less than 


Lemma 5.2. Let 0 <d < and let the grid boxes 
in D* have width w satisfying 


d 

W < -;=. 

" (Co + Ci)V2 


( 20 ) 


^ <Ci. 
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Then the following properties hold. 

1. The angle variation of Vh. over any gridbox is 
less than -0. 

2. Let p be a point on an edge of a gridbox, and 
let be the point on the boundary of the gridbox 
into which X^(p) is pointing, such that the line 
segment pq^ has direction Xa(p). The point q_^ 
is defined similarly. See Figure 12. Then Vh. is 
pointing rightward along pq# and leftward along 

pq-«- 

3. The function h is strictly increasing on the line 
segments from p to q^ and from p to q-#. 




q^ 

^ t 

Vh f 


V 

T / 

kj-' 

V 

.q-« 


Ij 


Figure 12: The orientation of Vh with respect to 
X^(p) does not change over a grid box. 

Proof. The first claim follows from Lemma 5.1, using 
the fact that the diameter of a grid box is wy/l. 

With regard to the second claim, the small an¬ 
gle variation condition implies that the orientation 
of {Vh(q),Xa(p)} does not change as q ranges over 
I. Since this orientation is positive for q = p, it is 
positive for all q G I. Similarly, the orientation of 
{Vh(q),X_.a(p)} is negative for all q G I. Therefore, 
the second claim also holds. 

At a point r of the line segment pq±« the direc¬ 
tional derivative of h in the direction of this line seg¬ 
ment is (Vh(r),X±^{p)), which is positive since the 
angle between Vh(r) and X±a(p), is less than •&, and 
^ < j. This proves the third part. □ 

Fencing in the separatrices For each isolating 
unstable separatrix interval J on the boundary of a 


saddle box we construct two polylines I_and 

Le(JJ) as follows. The initial points of these polylines 
are the endpoints of J, q_ and q+, where q_ comes 
before q+ in the counterclockwise orientation of the 
boundary of the saddle box. The polyline Le{J) is 
uniquely defined by requiring that its vertices q+ = 
po,pi,...,Pn lie on grid edges, with the property 
that 

1. The line segment pi_ipi, 0 < i < n, lies in 
a (closed) grid box of D*, and has direction 

X^ (pi-i). 

2. prt, the last vertex, lies on the boundary of D*. 

The polyline I_#{JI) is defined similarly, with the ob¬ 

vious changes: its initial vertex is q_, and each edge 
has direction equal to the value of the vector field X_^ 
at the initial point of this edge. The polylines L±^{JI) 
are called fences of the (unique) unstable separatrix 
of Vh intersecting J. 

It is not hard to see that that a grid box contains 
at most two consecutive edges of each of these poly¬ 
lines, but it is not obvious a priori that each box 
cannot contain more than two edges of each polyline 
in total. It follows from the next result that the in¬ 
tersection of a grid box with any of these polylines is 
connected, and, hence, that these polylines are finite. 

The following results states that, when walking 

along the polylines and 1_^ in the direction of 

increasing h-values, each grid box is passed at most 
once. 

Lemma 5.3. Let I be a box such that the angle 
variation of Vh over the surrounding box N (I) is 
at most Then the intersection of L^(JJ) and 

I (I_#(JI) and ts either empty or a connected 

polyline (consisting of one or two segments). 

Proof. Let q be a point at which La(J) leaves I, i.e., 
the segment of (J) ending at q lies inside I and the 
segment qr beginning at q lies outside I. Let p be 
a point on the boundary of I at which h attains its 
maximum value M. 
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Figure 13: The value of h. at the point where polyline 
Le (J] leaves the surrounding box N (I) is greater than 
the maximum value of h. on I. 


Case 1: is a vertex of I, incident to the edge of 
I containing q. 

See Figure 13, top row. Let I be the line through the 
edge of I containing q, let a be the angle between 
I and Vh(p), and let (3 be the angle between I and 
the segment of L^(JJ) with initial point q. The angles 
a and |3 are both positive, since p is a vertex of I. 
The angle between Vh.(p) and Vh(q) is at most 5, 
since the angle variation of Vh. over I is less than 
Therefore, |a — |3| < 2-0. 

Let ppo and ppi, with po and pi on the boundary 
of N (I) , be the line segments that make an angle of 
^ ^ with Vh.(p). Since the angle variation of Vh. 

over N (I) is at most the value of h. at any point of 
these line segments is at least M. We shall prove that 
the connected component of fl N{I) containing 
q intersects one of the line segments ppo and ppi. 

First assume a > 0. Then the line segment ppo lies 
in the grid box J[ containing segment qr of Lo(JJ). If 
r lies on an edge of JJ incident to p, then qr intersects 
ppo- So assume r lies on the edge of JJ contained in 


the boundary of N(I) (Figure 13, leftmost picture). 
Let s be the point of intersection of the line through 
ppo and the line through qr. This point lies on the 
same side of 1 as po and r, since Zpopq = ^ — a + 
0 < 5 and Zpqr = |3 < j. Furthermore, Zpsq = 
7t—(j — a + 0) — |3 > 0> j. Therefore, x lies 

inside N(I), in other words, qr intersects ppo also in 
this case. 

Now consider the case a < 0. Then po lies on 
the side of N(I) parallel to the line through p and q. 
Furthermore, |3 < a+20 < 30 < so L^(J[) ‘leaves’ 
N (I) at a point t on the side of N (I) perpendicular to 
the line through p and q. See Figure 13, rightmost 
picture. It follows that the part of L^{JI) between q 
and t intersects ppo- In particular, h.(t) > M. 

Case 2: p is not a vertex of I, incident to the edge 
of I containing q. Then either p is a vertex of I, not 
incident to the edge of I, containing q, as in Fig¬ 
ure 13, bottom-left picture, or p lies on the relative 
interior of an edge of I, as in Figure 13, bottom-right 
picture. 

In this case Vh.(p) is nearly vertical, as are the 
edges of Lo(JJ). Similarly, the line segments ppo and 
ppi are nearly horizontal, so Lo(J) intersects ppo- 
The details are similar to those of Case 1 of this proof. 

□ 


If the endpoints of the fences Lo(JJ) and I_o(J[) lie 

on the same connected component of the boundary 
of ID*, then these fences split ID* into two connected 
regions. See Figure 14. In this case, the region con¬ 
taining the separatrix interval JI in its boundary is 
called the funnel of JJ (with angle 0) denoted by 
F«{JI). Its boundary consists of JJ, the two fences 

Ls(JJ) and I_^(JJ), and a curve JJ* on the boundary of 

0D* connecting the endpoints of these fences. If the 
funnel is simply connected, it contains the part of the 
unstable separatrix through JJ lying inside D*, which 
enters the funnel through JJ and leaves it through JJ*. 
Note that JJ* is a curve either on the outer boundary 
of ID* or on a sink box. 
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Figure 14: Fences around a separatrix y. If the fences 
end in the same connected component of the bound¬ 
ary of mi*, then they enclose a funnel (top right pic¬ 
ture). If the funnel is simply connected, it isolates the 
separatrix from the source-, sink- and saddle-boxes 
(bottom picture). 


Similarly, each stable separatrix interval has two 
fences (for an angle ^). If the endpoints of these 
fences lie on the same connected component of 9D*, 
the enclosed region is again called a funnel for the 
stable separatrix interval. Our goal is to construct 
disjoint, simply connected funnels for the stable and 
unstable separatrix intervals. If these funnels are dis¬ 
joint, then they form, together with the sink boxes, 
source boxes and saddle boxes, a (fattened) Morse- 
Smale complex for Vh. 


It is intuitively clear that a funnel F^(J) is sim¬ 
ply connected if the length of the separatrix in¬ 
terval J, and the edge length w of the grid boxes 
are sufficiently small. The next subsection presents 
computable upper bounds on these quantities, guar¬ 
anteeing that the endpoints of two fences of a sepa¬ 
ratrix interval lie on the same boundary component. 
It is then easy to check whether the enclosed funnel 
is simply connected. 


5.2 Controlling the width of the fun¬ 
nel 

If the width of a funnel is sufficiently small, in a 
sense to be made more precise, it encloses a simply 
connected region in D*. The width of a funnel is, 
roughly speaking, the number of grid boxes between 
its bounding fences in the vertical direction, in re¬ 
gions where the fences are nearly horizontal, and in 
the horizontal direction, in regions where the fences 
are nearly vertical. To define the width of a funnel 
more precisely, we distinguish quasihorizontal and 
quasivertical parts of a funnel, and show that the 
width of a funnel does not increase substantially at 
transitions between these quasihorizontal and qua¬ 
sivertical parts. 


Quasihorizontal and quasivertical parts of a 
funnel A nonzero vector v = (vi,V 2 ) is called 
quasihorizontal if |v 2 | < 2|vi|, and quasivertical 
if |vi| < 2 |v 2|. Note that each nonzero vector is 
quasihorizontal, quasivertical, or both. Consider a 
subdivision of D* into boxes of equal width, where 
non-disjoint boxes share either an edge or a vertex. 
A horizontal e-strip is the union of a sequence of 
boxes where successive boxes share a vertical edge, 
such that the horizontal edge of the rectangle thus 
obtained has length at most e. A vertical e-strip is 
defined similarly. An e-box is a square box with edge 
length at most e which is the union of a number of 

boxes. Two polygonal curves L+ and I_form an e- 

funnel if there is a set "He of horizontal £-strips, a set 
Ve of vertical £-strips, and a set Se of £-boxes such 
that the following holds: 

1. The vertices of 1_and L+ lie on the edges of the 

grid-boxes; I_intersects a grid box in at most 

one vertex or in at most one edge; the same holds 
for L+; 

2. Both 1_and L+ lie in the union of the rectangles 

in 'He, Ve and Bt] 
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3. An edge of L± contained in a horizontal £-strip 

is quasi-vertical, and an edge contained in a ver¬ 
tical £-strip is quasihorizontal. Moreover, nei¬ 
ther L+ nor I_intersect the vertical sides of a 

horizontal £-strip, or the horizontal sides of a 
vertical £-strip. Each £-strip and each £-box is 
intersected by both polylines. 

4. Either I_or L+ intersects an £-box in exactly 

one of its edges, which is contained in a grid box 
at the corner of the £-box. This single edge is 
either quasivertical or quasihorizontal (but not 
both). If this edge is quasihorizontal (quasiver¬ 
tical), all edges of the other polyline inside the 
£-box are quasihorizontal (quasivertical) as well 
- and possibly also quasivertical (quasihorizon¬ 
tal). The other polygonal curve intersects the 
same edges of the £-box, each in exactly one 
point, and is disjoint from the other edges of 
the £-box. 


See also Figure 15. 

We determine ^ > 0 later, but for now we assume 


that 


^ < 


7t 

40 ■ 


( 21 ) 


We start with a simple observation. 


Lemma 5.4. Let L be a polyline with quasihori¬ 
zontal edges and with vertices on the edges of a 
grid with edge length w satisfying (20). If L lies 
in a vertical strip of width w, where each of the 
vertical lines bounding the strip contains one of 
its endpoints, then L intersects at most three grid 
boxes contained in this vertical strip. 

A similar property holds for a polyline with qua¬ 
sivertical edges intersecting a horizontal strip. 


Proof. We only prove the first part, in which L lies 
in a vertical strip and has quasihorizontal edges. The 
slope of the line segment connecting the endpoints of 
L does not exceed the maximum slope of any of the 
edges of L, so this slope is at most arctan2. Hence the 
projection of this line segment on any of the vertical 



Figure 15: A funnel formed by two polylines covered 
by two vertical £-strips, one £-box and three horizon¬ 
tal £-strips. Here £ is six times the width of a grid 
box. L+ intersects the £-box in a single edge, which 
is quasivertical but not quasihorizontal. All edges 

of I_inside the £-box are quasivertical as well (and 

some of them are also quasihorizontal). 


lines bounding the strip has length at most 2w, so it 
intersects at most three boxes. □ 

The next result shows that the width of a funnel 
does not grow substantially at a transition between 
a quasihorizontal and a quasivertical part. We take 
£ > 0 such that the angle variation of Vh. over a 
box with edge length £ is at most Again, by 
Lemma 5.1, this is guaranteed by taking 

£ < - - -( 22 ) 

“20(Co + C,)V2 

Lemma 5.5. Let J be an e-box intersected by both 

I_ 0 o-'n-d L+#, with an edge e which contains the 

initial vertex of both I_^ H J and L+# fl J. Assume 

that at least one of the polylines has an edge which 
is either quasihorizontal or quasivertical, but not 
both. Then both polylines intersect the boundary 
of J in exactly two points, and there is an edge e' 
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of I, adjacent to e, containing the terminal ver¬ 
tices of both I_^ n J and L+# fl J. See Figure 15. 

Proof. Assume that L_|_^ has an edge e+ which is 
quasivertical but not quasihorizontal. We first show 
that all edges of are quasivertical (and possibly 
quasihorizontal). The angle between e+ and the hor¬ 
izontal direction is at least arctan2, which is greater 
than 7 + Since the slope of e+ is the slope of 
the vector field X® at the initial vertex of e+, and 
the angle variation of over J is at most the 
slope of X® at any point of J is at least j + Since 
the slope of an edge of is the slope of X^ at the 
initial vertex of this edge, we conclude that all edges 
of are quasivertical. 

All edges of 1_^ are also quasivertical (and pos¬ 

sibly quasihorizontal). To see this, observe that the 

slope of an edge of 1_^ is the slope of X_# at the 

initial vertex of this edge, and, hence, the slope of 
X^ at this initial vertex, minus 2^. In other words, 

the slope of any edge of I _# is at least j + ^ — 2.d. 

Since ^ this slope is at least j. Therefore, all 

edges of I _^ are quasivertical. 

The polylines L+# and I_a do not intersect the 

edge of J opposite e, since then at least one of the 
edges of these polylines would have a slope less than 
j. Let e' be the edge containing the endpoint of 
La (J) nJJ. Then e' is adjacent to e. Given the bounds 
on the slope variation of the edges of the polylines, 
it is easy to see that 

(i) the endpoint of La(JJ) is the only point of this 
polyline on e'; 

(ii) the endpoint 1 _a H J also lies on e', and this is 

the only point of this polyline on e'; 

(iii) none of the polylines intersects the edge opposite 
e'. 

This concludes the proof of Lemma 5.5. □ 


bounds are conservative, they provide the tools for 
the construction of certified funnels for all separatri- 
ces. 

A gridbox is called quasihorizontal {quasiverti¬ 
cal) if it contains a point at which Vh. is quasihori¬ 
zontal (quasivertical). Again, a gridbox may be both 
quasihorizontal and quasivertical. 

An integral curve of Vh. in a quasihorizontal grid- 
box [xo,xi] X [r)o,yi] is the graph of a function 
X I— > ij(x), where y(x) is a solution of the differen¬ 
tial equation 

h'W = F(x,y(x)), (23) 

ij(xo) =yo, 


where LlXjy) = 




Here x ranges over the full 
Otherwise, the 


interval [xo,xi] if ijo < ij(x) <111 
range of x is restricted to a suitable maximal subin¬ 
terval [£,o,f.i], such that (f.o,y(f,o)) and (£.i,y(f,i)) 
are points on the boundary of the gridbox. Simi¬ 
larly, a trajectory of X^ in a quasihorizontal gridbox 
[xo, xi] X [ijo , y 1 ] is the graph of a function x i—> y (x), 
where y(x) is a solution of the differential equation 


dx 


= F^5(x,y), 


(24) 


with 


F«(x,y) = 


hx(x,y) sinff + hy(x,y) cosff 
hx(x,y) cosff — hy(x,y) sinff’ 


Similarly, a trajectory of X_^ is the graph of a 
function y x(y), where x(y) is a solution of the 
differential equation 


with 


dx 


GD(x,y), 


Growth of the width of quasihorizontal and 
quasivertical funnel parts The width of the fun¬ 
nel may grow exponentially in the number of grid 
boxes it is traversing. The next result gives an upper 
bound for the growth of this width. Even though the 


g 1 _ hx(x,y) cosff-hy(x,y) sinff 

^ F#(x,y) hx(x,y) sinff + hy(x,y) COS-0 

Here y ranges over the full interval [yo,yi] if xq < 
x(y) < xi, or a suitable maximal subinterval other¬ 
wise. 
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The union of all quasihorizontal gridboxes in D* is for |^| < ^qy. Finally, let the constants Co, Ci and 
denoted by ID>*qh, and the union of all quasivertical be defined by 

gridboxes by D* Co = 2max(Cqh + AqhBqi,, Cqy + AqyBqy) (29) 

Even though the width of a funnel may grow ex- ^ ^ 

ponentially in the number of grid boxes it traverses, Ci = max(-,-) (30) 

this growth is controlled. To this end, we introduce 2MqhKqii 2MqyKqy 

several computable constants that only depend on = niax(5qh,5qv)- (31) 

the function h. and (the size of) its domain ID . Let result provides an upper bound for the 

Aqh, Aqy, Bqii, Bqy, Cqh and Cqv be positive con- growth of the funnel width along a quasihorizontal 

stants such that bounding polylines. We assume that the 

IT-, M ^ A IT-, II ^ A funnel runs from left to right, so its initial points are 

max |r(x,yj| < Aqh, max lC3lx,yj| < Aqy, . . . 

(x,y)eB*qh (x,y)GD*qv on the line with smallest x-coordmate. If the funnel 

,T- -AT- runs from right to left, a similar result is obtained, 

of oG, Tj 

5.6. Let : [a,b] -> [c, d] be the 

piecewise linear functions the graphs of which are 


max |F(x,y)| < Aqh, 
(x,y )GB*qh 

max 

(x,y)GD- 

max 

(x,y)GD*qh 

l-^Ky)! < Bqh, 

max 

(x,y)GD" 

max 

(x,y )GB*qh 

91” 

l^(x,y)| < Cqh, 

max 

(x,y)GD* 

Note that 




c,y)| < parts of the polylines and L_ 


Fo(x,i))-F(x,'y) = 


(Bxlx,^)^ -hh,y{x,y)2) sin^ 


for a grid with edge length w, respectively. Let A 
be an upper bound for the distance of the initial 
points of these polylines, i.e., 


h.x(x,'y)^ COS-0 — hxlXjij) by (Xjy) sin0 


lyB,w(Q) -y-B.-wla)! < A. 


Let Mqii be a dyadic number such that 


By(x,l)) 


max -- 

(x,y)GD*,h hxlXj-y) 


^ '-qh • 


Then the width of the fence, bounded by and 
I_ is bounded: 

gCqh(x a) - 1 

(25) |ij^,Tv(x)-y-.9,Tv(x)| < Ae‘“‘>*'^’‘^“’-|-(cow-|-ci sin0)--- 

Gqh 


1 Proof. Let y±^(x) be the exact solution of the ro- 

Take 0qh € (0, jTt) such that tan0qh < j • Fi- tated system with initial condition ij±^(a). In par- 

, 2 , ticular, |'y#(a) — ■y_^(a)| < A. Then (27) implies 


nally, let be a constant such that 


h.x(x,v)^ -hhy(x,'q)^ 


-(x)-F(x,'y±B(x))| = |F±B(x,-y±B(x)-F(x,y±B(x))| < Mqhsinb 


(x,y)GD*qh| hxlXj-y)^ COS-0qh 


^ - V • 


Therefore, according to the Fundamental Inequal¬ 
ity [16, Theorem 4.4.1]—See also A—we have 


Taking Mqh = —itT) have, for \d\ < -^qh: 


|yi)(x)-y_^)(x)| < ). 

Cqh 


rnax |F#(x,y) F{x,y)| < Mqh sin0. (27) The interval [a, b] is subdivided into a finite number 

{x,y)GB*qh 

of subintervals of length at most w, where the end- 


Similarly, there are (computable) constants Mqy and points correspond to the x-coordinates of the break- 
^qv such that points of the fences and I_Let ya,w be the Eu¬ 

ler approximation to the ordinary differential equa- 
(x,™eD*qy G(x,y)| < Mqv sinff, (28) tion (24). Its graph is (a quasihorizontal) part of the 
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fence L^. Theorem 4.5.2 in [16]—See also A—gives 
the following explicit bound for the error in Euler’s 
method: 

^qh 

(33) 

We get a similar inequality for the Euler approxima¬ 
tion of Combining (32) and (33), and 

using (29) and (30), yields 

IU«,w(x) -■y_^),w(x)| < 

2(w(Bqh + AqhCqh) + Mqh sin-0) 

= A + (cqw + Cl sin-O) 

□ 

A similar result holds for quasivertical trajectories. 
Next we need to control the increase of the funnel 
width upon transition from a quasihorizontal to a 
quasivertical part its bounding polylines (or from a 
quasivertical to a quasihorizontal part). 

Transitions: bounded increase of funnel width 
Transition from a quasihorizontal to a quasivertical, 
or from a quasivertical to a quasihorizontal part of 
the funnel takes place at an e-box. If the width of 
the funnel at the ‘entry’ of the box is less than the 
width w of a grid box, then the width may increase, 
but it will not be greater than 2w at the exit. This 
is made more precise by the following result. 

Lemma 5.7. Let I be a e-box as in Lemma 5.5, 
where, moreover, the initial points p and q o/ fl 

J and I_ 0 n J, respectively, are on the boundary 

of the gridboxes containing the vertices of edge e 
of I. If the distance between p and q is at least 
w, then the distance between the terminal points 

p and ■q o/ fl JJ and I_^ fl J, respectively, is 

less than the distance of p and q. If the distance 
between p and q is less than w, then the distance 
o/p and ■q is at most 2w. 


Proof. Assume that the first edge of is quasiver¬ 
tical, but not quasihorizontal. Edge e of JJ is then 
vertical. Assume that this polyline consists of a sin¬ 
gle edge, namely the line segment pp. Let |3+ be the 




-a] _ I 



d 

V q 



J?igilBe-46;_'^he distance between the two polylines 
upor^gtry and exit of a box. Left: If the distance d 
between the initial points p and q of the polylines is 
less than the edge-length w of a grid box, then the 
distance d between the terminal points p and q' is less 
than 2w. Right: Otherwise, the distance d between 
the terminal points is less than d. 


angle between pp and edge e, then arctan j < |3+ < 
^ -|- arctan 2 . Let |3_ be the angle between the line 
segment q'q and edge e, then |3_ is inbetween the 

smallest and largest slope of any edge of I_a- Since 

the angle variation of X over J is less than the 
angle |3- is greater than |3+ — Let a be the dis¬ 
tance of p to the nearest vertex of e, then a < w. If 
d > w, the distance d between p and q' satisfies 

d = (d -|- a) tan |3- — atan |3+ 

< d tan(|3+ -h ^) -f a (tan(|3+ + ^) — tan |3+) 

< |d-h la 

<d, 


since a < w < d. Here we used tan |3+ < j to get 


tan(|3+-|-^) 


tan |3+ -F tan ^ ^ 1 + tan ^ 3 
1 — tan |3+ tan ~ I ~ j tan ^ ~ 


Since arctan j — ^ < arctan j — ^ < |3+ < arctan j, 
a short computation shows that tan(|3+ + 2^") - 
tan|3+ < 1. 
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If d < w, then q lies in the same gridbox as p, or in 
a gridbox adjacent to it. Then it is easy to see that p 
lies in the same grid box as p, and q' also lies in this 
box, or in a box adjacent to it. Therefore, d < 2w in 
this case. 

If I_^ consists of a single edge, then the argument is 

similar. □ 


does not exceed A^-i, cf Lemma 5.7. Assume that 
the n-th part of the funnel is quasihorizontal, then 
Lemma 5.6 implies that the width of this part at a 
point with horizontal coordinate x is at most 

r f 1 eCqh(x-a)_i 

An -1 + (cow + Cl sin^)---, 

^qh 


T r- r 1 r .7 • j j .1 r 11 • ,, SO in particular, since D < Cnh < C and 0 < x — a < 

Lemmas 5.6 and 5.7 provide the following result ^ — — — 

on the upper bound on the funnel width of a separa- 

trix with M transitions between quasihorizontal and a ^ a ct , Cow + Cisin^, 

, An < An-l e H -- jC — IJ. 

quasivertical parts. U 


Corollary 5.8. Let T he the (computable) edge 
length of a bounding square of the domain D of 
the function \\), and let M he the total num¬ 
ber of quasihorizontal and quasivertical parts of 
the polylines bounding a separatrix funnel. Let 
C = max(Cqh,Cqv) and let D = min(Cqh, Cqv). 
Then the width of the funnel does not exceed 

gCMT 

(Cld + C2W] ^ , 

Co 

provided d < do where ci “ ^ "q > with Co and ci 
given by (29) and (30), respectively. 

In particular, this width is at most £ if 

Cid + C2W<(34) 

Proof. Let D* C [a, b] x [c, d], then T < max(b — 
a, d — c). There are M — 1 transitions from quasi¬ 
horizontal to quasivertical parts of the funnel, each 
occurring at an e-box. Let Ao be the width of the ini¬ 
tial separatrix interval, and let Ai,..., Am-i be the 
width of the funnel at the entry of the corresponding 
boxes, in other words, A^ is the width at the end of 
the k-th part of the funnel. Using induction, we will 
prove that, for k = 1,..., M: 


A ^ T kCT Cow -|- Cl sind , 

Ak < 2w e H-—-(e — 1, 




So assume (35) holds for k = n — 1. If An-l > 
w, the initial width of the n-th part of the funnel 


Therefore, (35) holds for k = n. If An-i < w, then 
the initial width of the n-th part of the funnel is 
at most 2w, cf Lemma 5.7. Therefore, Lemma 5.6 
implies that the width of this part at a point with 
horizontal coordinate x is at most 


-h (cqw- h Cl sind) 


so in particular 


jCqh(x-a) _ 1 


c, 


qh 


< iwe'C'’ + - I), 

Therefore, for n = M, we have 


Am < (ci sind + C 2 w)- 


,CMT 


- 1 


D 


< (Cid -h C2W) 


,CMT 


D 


which proves the corollary. 


□ 


Remark 5.9. The computable constants 
do, Cl, C 2 , C and D depend only on D* and h. 

□ 


In the next section, we assemble the bits and pieces 
into a certified algorithm for the construction of the 
MS-complex, and show how the upper bounds on the 
funnel width are used to prove that this algorithm 
terminates. 
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5.3 Construction of the MS-complex 

The Algorithm The construction of the MS- 
complex of h is a rather straightforward application 
of the preceding results. It uses a parameter M, the 
(a priori unknown) number of transitions (at e-boxes) 
between quasihorizontal and quasivertical parts of a 
funnel. Let T be the edge length of a bounding square 
of the domain D of h. Then the algorithm performs 
the following steps. 

Step 1. Construct certified isolating boxes 
,..., for the singularities of Vh (cf Sec¬ 
tion 3). 

Step 2. Let ID* be the closure of D \ {B( U • • • U 
B(^). Compute the constants -Do, Ci, Ci, C and 
D, which depend only on h and ID* . Set e to the 
minimum of the width of the source-, sink- and 
saddleboxes. 

Step 3. Let ^ and w be such that w < 

■0 < min(^,'0o), and ci-Q + ciw < e D 
(cf Corollary 5.8). Subdivide D* until all grid- 
boxes have maximum width w. For each sad¬ 
dle box, compute four separatrix intervals on its 
boundary, of width at most w. 

Step 4. For each stable and unstable separatrix in¬ 
terval do the following. Start the computation 
of a funnel for a separatrix by setting M to a 
small number Mo (say 4). Compute the fences 

I_ 0 and Lo, keeping track of the width of the 

enclosed funnel under construction and of the 
number m of transitions between quasihorizon¬ 
tal and quasivertical parts of this funnel. 

If the width of the funnel exceeds e or the 
number of transitions m exceeds M, then 
abort the computation of the current funnel, dis¬ 
card all funnels constructed so far, set M to twice 
its current value and goto Step 3. 

If the funnel intersects an already constructed 
funnel, or a source- or sinkbox on which it 
does not terminate (i.e., if only one of its fences 


intersects this box), then set £ to half its current 
value, discard all funnels constructed so far, and 
goto Step 3. 

If the funnel intersects a saddlebox B(, then 
decrease the size of B( by a factor of two via 
subdivision, discard all funnels constructed so 
far, set e to half its current value, and goto Step 
2. (Note that D* gets larger, so the constants in 
Step 2 have to be recomputed.) 

Otherwise, the fences end on the same compo¬ 
nent of the boundary of 0D*. The enclosed fun¬ 
nel is simply connected, and does not intersect 
any of the funnels constructed so far. Add the 
funnel to the output, and reset M to Mo (and re¬ 
peat until all separatrices have been processed). 


5.4 Termination 

Since the gradient field Vh is a 2D Morse-Smale sys¬ 
tem, its separatrices are disjoint. Their intersections 
with ID* are compact, and have positive distance (al¬ 
though this distance is not known a priori). In the 
main loop of the algorithm, the maximal funnel width 
£ is bisected if funnels intersect, and saddleboxes in¬ 
tersected by the funnel are subdivided, so after a fi¬ 
nite number of iterations of the main loop its value 
is less than half the minimum distance between any 
pair of distinct separatrices, and funnels stay clear 
from saddleboxes (apart from the one containing the 
a- or co-limit of the enclosed separatrix). 

Separatrices that intersect 01D do so transversally, 
cf Remark 2.1. Therefore, after a finite number of 
subdivision steps, both fences around such separa¬ 
trices will intersect 0D transversally. Hence, eventu¬ 
ally all funnels become disjoint, at which point the 
algorithm terminates after returning a topologically 
correct MS-complex for Vh. 
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6 Implementation and experi¬ 
mental results 

The algorithm has been implemented using the Boost 
library [1] for lA. All experiments have been per¬ 
formed on a 3GHz Intel Pentium 4 machine under 
Linux with 1 GB RAM using the g++ compiler, ver¬ 
sion 3.3.5. Figures l(a)-(b) and 17(a)-(b) depict the 
output of our algorithm, for several Morse-Smale 
functions. In our implementation the parameter 
used in the construction of separatrix-funnels, is 
which is larger than the theoretical bound given by 
Corollary 5.8. The algorithm halves this angle sev¬ 
eral times, depending on the input, until the funnels 
are simply connected, mutually disjoint, and connect 
a saddle-box to a source-box (for stable separatrices) 
or sink-box (for unstable separatrices), in which case 
a topologically correct MS-complex has been com¬ 
puted. Each of the funnels with deep black bound¬ 
aries contains an unstable separatrix, whereas a fun¬ 
nel with light black boundaries contains a stable sep¬ 
aratrix. The CPU-time for computing a MS-system 
increases with the number of critical points and the 
complexity of the vector field, as indicated in the 
captions of the figures. 

7 Conclusion 

The outcome of our research is two-fold. Firstly, 
we compute the topologically correct MS-complex of 
a Morse-Smale system. The correct saddle-sink or 
saddle-source connectivity can also be represented as 
a graph, which is of special interest from different 
application point of view. On the other hand, de¬ 
pending on a user-specified width of funnel one can 
compute a geometrically close approximation of the 
MS-complex. We give the proof of convergence of our 
algorithms. Although the complexity of the given al¬ 
gorithm depends on the input function and the com¬ 
plexity of the interval arithmetic library used in the 
algorithm. As we discussed some of the separatri- 



(a) Contour plot (left) and MS-Complex (right) of 

h(x,y) = — lOx^ -I- x"* + — y'^ + X + xy^, on the box 

[—4,3.5] X [—4,3.5]. CPU-time: 11 seconds. 


(b) Contour plot (left) and MS-complex (right) of a prod¬ 
uct of seven linear functions, on the box [—7,7] X [—7,7]. 
CPU-time: 11.5 minutes. 

Figure 17: Contour plots of MS-functions and their 
Morse-Smale complexes. 

ces inside a bounding box B may have discontinuous 
components. The algorithm we propose here is able 
to compute only the part of the separatrices which 
are connected to the corresponding saddle. Therefore 
one open question is how to compute all the compo¬ 
nents of separatrices inside a bounding box. 
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Error in Euler’s method. Error bounds for ap¬ 
proximate solutions of ordinary differential equation 
play a crucial role in the construction of certified fun¬ 
nels for separatrices. We quote the relevant parts of 
the book [16]. 

Fundamental Inequality [16, Theorem 4.4.1]. 
Consider the differential equation 


dx 




on a box B = [a, b] x [c, d], and let C be a constant 
such that 


max 

(x,y)eB 




c. 


7/yi(x) and yzix] are two approximate piecewise 
differentiable solutions satisfying 


|y[(x) -F(x,'yi[x))| < £i, 

|y2W -F(x,'y2W)| < £2 


for all X G [a, b] at which yi (x) and yilx) are dif¬ 
ferentiable, and if, for some xo G [a, b] 
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lyi (xo) -yilxo)! < 5, 
then, for all x G [a, b] 

pC|x-xo| _ 1 

lyi(x) -y2(x)| < -he- - -, 

where e = £i -h £ 2 - 

The well-known Euler method for constructing ap¬ 
proximate solutions to ordinary differential equations 
is also useful for the construction of certified strips. 
It proceeds as follows. For a given initial position 
(xo,yo), define the sequence of points (xn,yn) by 

Xn =Xn-i -hri =xo -hnq 
hn =hn-l +riF(Xn-l,yn-l), 
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as long as (xTi,yri) S B. Then the Euler approxi¬ 
mate solution through (xo,i)o] with step r| is 

the piecewise linear function the graph of which joins 
the points (xr^ijn), so 

Uh(x) =-yn + (x-Xn)F(XTx,yn) for X £ [Xn,Xn+l]. 

The following result states that the Euler approxi¬ 
mate solution converges to the actual solution as the 
step tends to zero, and gives a bound for the error. 

Error in Euler’s method [16, Theorem 4.5.2]. 
Consider the differential equation 

on a box B = [a, b] x [c, d], where T is a -function 
on B. Let the constants A, B and C satisfy 

9P 

max |F[x,y)| < A, max | —(xj-y)] < B, max 

(x,y)eB OX (x,y)el 

The deviation of the Euler approximate solution 
yri with step y from a solution y of the differential 
equation with lyla] — ^^ 1 ( 0 )! < A satisfies 



x = bo x = bi x = b 2 

Figure 18: Zooming in on the saddle point by subdi¬ 
vision. 

2. the saddle point p is contained in box In, for all 
n. 

See Figure 18. Let s be the x-coordinate of the sad¬ 
dle point p, and let bn be the x-coordinate of the 
rightmost vertical boundary edge of N(In]. Let Wn 
l^^the width of In, and let Cn be the x-coordinate of 
|i|s-^8nt)^.<r£en bn = Cn + fwn, and Wn+l = jWn. 
since then b = bo > bi > ..., since 

bn+l = Cn+l+fwn+1 < Cn + {wn + fwn = bn-yWn. 

Since |cn — s| < -jWn, we get 


lyri(x) -y(x)| < “'+y(B + AC) 


,C|x-a|-l 

c 


Wn < bn — S < 2 Wn. 


(36) 


for all X G [a, b]. 

The preceding result also holds if, as in the current 
chapter, y is not the exact step, but an upper bound 
for a possibly varying step. 

B Narrowing separatrix inter¬ 
vals 

We first sketch the algorithm for narrowing the sepa¬ 
ratrix intervals. To this end we subdivide the box I, 
and hence the box N(I), yielding a nested sequence 
of boxes I = lo A Ii A ..., with surrounding boxes 
N(I) = N(Io] A N(Ii) A ..., such that 

1. width(In+i ] = 2 width(In) 


Consider the forward integral curves of the vector 
field Vh through the points of intersection qn of the 
line X = bn and the boundary curves F^p. See Fig¬ 
ure 19. These curves intersect the rightmost edge of 
N(I) in two points bounding an interval Jn on this 
edge. Arbitrarily good approximations of these in¬ 
tegral curves are obtained as follows. Let ^n be (an 
upper bound on) the maximum angle variation of Vb 
over any of the boxes of the n-th subdivision of N (I). 
Since b is , the angle variation is a Lipschitz func¬ 
tion, so limn-joo^n = 0- In particular, the rotated 
vector fields X±#. converge to Vb. We construct an 
upper fence with angle for the upper integral 
curve, and a lower fence with angle — jbn for the 
lower integral curve. See also Section 5.1 for the con¬ 
struction of a fence. These fences are disjoint, since 
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the angle variation of Vh. over a grid box is less than 

^n. 

Since limn->oo bn = s, cf (36), the points con¬ 
verge to the saddle point. Therefore, the intervals 
Jo D Ji D . •contained in the intersection of the 
unstable wedge Cp and the rightmost vertical edge 
of N (I), converge to the intersection of the unstable 
separatrix and the rightmost vertical edge of N (I). 



A proof of Lemma 4.4 can be given along the lines 
of [9, page 330ff] or [8, Chapter 3.6]. Rather than 
giving a complete proof we give an example illustrat¬ 
ing the main ideas. Consider the function h{x,y) = 
-h 2^5 with As < 0 < Au. The gradient 
vector field is given by Vh{x,y) = (A^x, Ajy)^. Ob¬ 
viously, the origin is a saddle point, with positive 
eigenvalue A^^ and negative eigenvalue Aj, and eigen¬ 
vectors (1,0)^ and (0,1)^, respectively. The unsta¬ 
ble cone of this saddle point is bounded by the curves 
defined implicitly by r))^p(x,y) = 0, where 

^±p(x,y) = det(V’^,Xp(x,y)) 

1 (Au cos |3) X =F (As sin |3) y 
0 (Au sin (3) X ± (As cos (3) y 

= (Au sin (3) X ± (As cos (3) y. 


Therefore, the equation of is y = ±ax, with 
a = —(tan|3)^ > 0. Let the right vertical edge of 
I be on the line x = b, b > 0, and consider a point 
q+ = (f,, af,), with 0 < £, < b, on the boundary 
curve Fp of the unstable cone. The integral curve of 
Vh. through q+ satisfies the differential equation 



y 

X 


with initial condition y (£,) 
Therefore, 


a£., where A = — < 0. 

Au 


y(x) = a£.(^)'^. 


The integral curve through q+ intersects the right¬ 
most edge of I in the point (b, 5 (£,)), where 


6(£.) = 


ab^ 


Similarly, the integral curve through q^ = 
(£,, —a£,) G F^p intersects the rightmost edge of I 
in the point (b, —6(f,)). Now let £, range over the 
sequence bo,bi,_ Then the interval Jn has end¬ 

points (brn±6(bn)), so its width is 25(bn). In view 
of (36), with s = 0, we have 

ab'^ . , WTT ^ -I 

^ A-1 ^ width(Jn) < 2—- 

Wn [2-Wnr 


In other words, with K = 1 — A > 1 and c = ab'^Wg , 
c < width(Jn) < 2’^ c 


Hence, 

width(Jn+ 2 ) < ^width(Jn). 

Since K > 1, after two subdivision steps the size of 
the separatrix interval reduces by more than a fac¬ 
tor two. Hence interval arithmetic provides an arbi¬ 
trarily good approximation of the intersection of the 
unstable separatrix and the boundary of the saddle 
box. A similar observation holds for the intersection 
of the other separatrices and the boundary of their 
saddle box. 
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